#define PY_ARRAY_UNIQUE_SYMBOL pbcvt_ARRAY_API

#include <boost/python.hpp>
#include <pyboostcvconverter/pyboostcvconverter.hpp>


#include <opencv/cv.h>
#include <boost/filesystem/operations.hpp>
#include <gtest/gtest.h>
#include <glog/logging.h>
#include <opencv2/opencv.hpp>
#include <stdbool.h>

#define PI 3.1416


namespace pbcvt {

    using namespace boost::python;
    using cv::Mat;
    using std::vector;

/**
 * Example function. Basic inner matrix product using explicit matrix conversion.
 * @param left left-hand matrix operand (NdArray required)
 * @param right right-hand matrix operand (NdArray required)
 * @return an NdArray representing the dot-product of the left and right operands
 */
    PyObject *dot(PyObject *left, PyObject *right) {

        cv::Mat leftMat, rightMat;
        leftMat = pbcvt::fromNDArrayToMat(left);
        rightMat = pbcvt::fromNDArrayToMat(right);
        auto c1 = leftMat.cols, r2 = rightMat.rows;
        // Check that the 2-D matrices can be legally multiplied.
        if (c1 != r2) {
            PyErr_SetString(PyExc_TypeError,
                            "Incompatible sizes for matrix multiplication.");
            throw_error_already_set();
        }
        cv::Mat result = leftMat * rightMat;
        PyObject *ret = pbcvt::fromMatToNDArray(result);
        return ret;
    }

//This example uses Mat directly, but we won't need to worry about the conversion
/**
 * Example function. Basic inner matrix product using implicit matrix conversion.
 * @param leftMat left-hand matrix operand
 * @param rightMat right-hand matrix operand
 * @return an NdArray representing the dot-product of the left and right operands
 */
    cv::Mat dot2(cv::Mat leftMat, cv::Mat rightMat) {
        auto c1 = leftMat.cols, r2 = rightMat.rows;
        if (c1 != r2) {
            PyErr_SetString(PyExc_TypeError,
                            "Incompatible sizes for matrix multiplication.");
            throw_error_already_set();
        }
        cv::Mat result = leftMat * rightMat;

        return result;
    }


/*///////////////////////////////////////////////////////////////////////////////
 *  Vessel extract related functions
 */

    Mat find_connected(Mat img, int area = 1000);

    Mat getLineMask(int theta, int masksize);

    Mat rotatex(Mat mask);

    Mat getBaseMask(int theta, int masksize);

    Mat drawLine(cv::Point pt0, cv::Point pt1, Mat orgmask);

    Mat getLinePts(cv::Point pt0, cv::Point pt1);

    Mat grow_boundary(Mat img, Mat mask, int erosionsize = 21, int iterations = 20);

    Mat global_normalize(Mat img, Mat mask);

    Mat get_lineresponse(Mat img, int linescale, int w = 15);

    Mat create_mask(Mat _img, double thresh = 0.001);

//    Mat vessel_seg(Mat img, int index = 0);

    Mat myIm2Double(Mat in);

/*
 *  This function return CV_8UC1 Data
 */
    Mat find_connected(Mat img, int area) {
        Mat label, stats, centroids;
        connectedComponentsWithStats(img, label, stats, centroids, 4, CV_32S);
        Mat oimg(img.size[0], img.size[1], CV_8UC1);
        oimg = 0;
        for (int i = 1; i < stats.size[0]; i++) {
            if (stats.at<int32_t>(i, CC_STAT_AREA) < area)
                continue;
            for (int x = stats.at<int32_t>(i, CC_STAT_TOP);
                 x < stats.at<int32_t>(i, CC_STAT_TOP) + stats.at<int32_t>(i, CC_STAT_HEIGHT);
                 x++) {
                for (int y = stats.at<int32_t>(i, CC_STAT_LEFT);
                     y < stats.at<int32_t>(i, CC_STAT_LEFT) + stats.at<int32_t>(i, CC_STAT_WIDTH);
                     y++) {
                    if (label.at<uint32_t>(x, y) == i) {
                        oimg.at<uchar>(x, y) = 255;
                    }
                }
            }
//        printf("handled %d, area %d\n", i, stats.at<int32_t>(i,CC_STAT_AREA));
        }
        return oimg;
    }

    Mat find_areaRange(Mat img, int minarea, int maxarea) {
        Mat label, stats, centroids;
        connectedComponentsWithStats(img, label, stats, centroids, 8, CV_32S);
        Mat oimg(img.size[0], img.size[1], CV_8UC1);
        oimg = 0;
        for (int i = 1; i < stats.size[0]; i++) {
            if (stats.at<int32_t>(i, CC_STAT_AREA) < minarea ||
                stats.at<int32_t>(i, CC_STAT_AREA) > maxarea)
                continue;
            for (int x = stats.at<int32_t>(i, CC_STAT_TOP);
                 x < stats.at<int32_t>(i, CC_STAT_TOP) + stats.at<int32_t>(i, CC_STAT_HEIGHT);
                 x++) {
                for (int y = stats.at<int32_t>(i, CC_STAT_LEFT);
                     y < stats.at<int32_t>(i, CC_STAT_LEFT) + stats.at<int32_t>(i, CC_STAT_WIDTH);
                     y++) {
                    if (label.at<uint32_t>(x, y) == i) {
                        oimg.at<uchar>(x, y) = 255;
                    }
                }
            }
//        printf("handled %d, area %d\n", i, stats.at<int32_t>(i,CC_STAT_AREA));
        }
        return oimg;
    }

    Mat myIm2Int(Mat in) {
        Mat result;
        if (in.depth() == CV_8U) {
            return in;
        }
        CHECK(in.depth() == CV_32F || in.depth() == CV_64F);
        in.convertTo(result, CV_8UC(in.channels()), 255.0);
        return result;
    }

    Mat vessel_seg(Mat img, int index) {
        img = myIm2Int(img);
        Mat mask_img = create_mask(img);
        vector<Mat> image_channels;
        split(img, image_channels);
        Mat gChannel = image_channels[1];
        Mat mask_gChannel = mask_img;
        gChannel = 255 - gChannel;    /* Í¼Ïñ»Ò¶È·´×ª£¬Ê¹Ñª¹Ü»Ò¶È¸ßÓÚ±³¾° */
        Mat dilated = grow_boundary(gChannel, mask_gChannel);
        Mat dilated_double = Mat::zeros(mask_img.rows, mask_img.cols, CV_32FC1);
        Mat dilated_norm = Mat::zeros(mask_img.rows, mask_img.cols, CV_32FC1);
        dilated.convertTo(dilated_double, CV_32FC1, 1.0 / 255.0);
        dilated_norm = global_normalize(dilated_double, mask_gChannel);
        Mat line_response = Mat::zeros(mask_img.rows, mask_img.cols, CV_32FC1);
        Mat line_response_norm = Mat::zeros(mask_img.rows, mask_img.cols, CV_32FC1);

        int scales = 0;
        int w = 21;
        int step = 2;

#pragma omp parallel for shared(dilated_norm) private(line_response,line_response_norm)
        for (int line = 1; line < 21; line += step) {
            line_response = get_lineresponse(dilated_double, line, w);
            line_response_norm = global_normalize(line_response, mask_gChannel);
            dilated_norm = dilated_norm + line_response_norm;
            scales++;
        }
        Mat avg_response = dilated_norm / (1 + scales);

        Mat vessels = Mat::zeros(img.rows, img.cols, CV_8UC1);
        double thresh = 0.48;
        for (int i = 0; i < img.rows; i++)
            for (int j = 0; j < img.cols; j++)
                if (avg_response.at<float>(i, j) > thresh) {
                    vessels.at<uchar>(i, j) = 255;
                }
        vessels = find_connected(vessels, 5000);
        return vessels;
    }

    Mat get_lineresponse(Mat img, int linescale, int w) {
        Mat avg_kernel(w, w, CV_32FC1, Scalar(1));
        avg_kernel = avg_kernel / (w * w);
        Mat avg_response;
        Mat linemask(linescale, linescale, CV_8UC1, Scalar(0));
        Mat linemask_double(linescale, linescale, CV_32FC1, Scalar(0));
        Mat line_strength;
        Mat line_response = (-200.0) * Mat::ones(img.rows, img.cols, CV_32FC1);
        filter2D(img, avg_response, img.depth(), avg_kernel);
        for (int theta = 0; theta < 180; theta += 15) {
            linemask = getLineMask(theta, linescale) * 1.0;
            linemask.convertTo(linemask_double, CV_32FC1, 1.0 / 255.0);
            linemask_double = linemask_double / countNonZero(linemask_double);
            filter2D(img, line_strength, img.depth(), linemask_double);
            line_strength = line_strength - avg_response;
            for (int i = 0; i < img.rows; i++)
                for (int j = 0; j < img.cols; j++) {
                    if (line_strength.at<float>(i, j) > line_response.at<float>(i, j))
                        line_response.at<float>(i, j) = line_strength.at<float>(i, j);
                }
        }
        return line_response;
    }

    Mat grow_boundary(Mat img, Mat mask_img, int erosionsize, int iterations) {
        Mat dilated = Mat::zeros(mask_img.rows, mask_img.cols, CV_8UC1);
        for (int i = 0; i < mask_img.rows; i++)
            for (int j = 0; j < mask_img.cols; j++) {
                if (i == 0 || i == mask_img.rows - 1)
                    mask_img.at<uchar>(i, j) = 0;
                if (j == 0 || j == mask_img.cols - 1)
                    mask_img.at<uchar>(i, j) = 0;
            }
        Mat erode_element = getStructuringElement(MORPH_ELLIPSE, Size(erosionsize, erosionsize));
        Mat erode_mask = mask_img.clone();
        erode(mask_img, erode_mask, erode_element);
        dilated = img.mul(erode_mask / 255);
        Mat oldmask = erode_mask.clone();
        int filter_rows[9] = {-1, 0, 1, -1, 0, 1, -1, 0, 1};
        int filter_cols[9] = {-1, -1, -1, 0, 0, 0, 1, 1, 1};
        Mat newmask = Mat::zeros(mask_img.rows, mask_img.cols, CV_8UC1);
        Mat outerborder = Mat::zeros(mask_img.rows, mask_img.cols, CV_8UC1);
        int pixelrow, pixelcol;
        int total, count;

        Mat diamond(3, 3, CV_8U, Scalar(1));
        diamond.at<uchar>(0, 0) = 0;
        diamond.at<uchar>(2, 0) = 0;
        diamond.at<uchar>(0, 2) = 0;
        diamond.at<uchar>(2, 2) = 0;

        for (int i = 0; i < iterations; i++) {
            outerborder = Mat::zeros(mask_img.rows, mask_img.cols, CV_8UC1);
            dilate(oldmask, newmask, diamond);
            outerborder = newmask - oldmask;
            for (int j = 0; j < outerborder.rows; j++)
                for (int k = 0; k < outerborder.cols; k++)
                    if (outerborder.at<uchar>(j, k)) {
                        total = 0;
                        count = 0;
                        for (int h = 0; h < 9; h++) {
                            pixelrow = j + filter_rows[h];
                            pixelcol = k + filter_cols[h];
                            if (pixelrow < outerborder.rows && pixelrow >= 0 \
 && pixelcol < outerborder.cols && pixelcol >= 0 \
 && oldmask.at<uchar>(pixelrow, pixelcol)) {
                                total = total + dilated.at<uchar>(pixelrow, pixelcol);
                                count++;
                            }
                        }
                        dilated.at<uchar>(j, k) = (uchar) (total / count);
                    }
            oldmask = newmask.clone();
        }
        return dilated;
    }

    Mat global_normalize(Mat img, Mat mask) {
        Mat norm_img(img.rows, img.cols, CV_32FC1, Scalar(0));
        Scalar mean, stdDev;
        int used_pixels = countNonZero(mask);
        double diff;
        meanStdDev(img, mean, stdDev, mask);
        stdDev.val[0] = stdDev.val[0] * std::sqrt(used_pixels / (used_pixels - 1.0));
        for (int i = 0; i < norm_img.rows; i++)
            for (int j = 0; j < norm_img.cols; j++)
                if (mask.at<uchar>(i, j)) {
                    diff = (img.at<float>(i, j) - mean.val[0]);
                    if (diff > 0) {
                        norm_img.at<float>(i, j) = diff / stdDev.val[0];
                    }

                }
        return norm_img;
    }


    Mat getLinePts(Point pt0, Point pt1) {
        int x0 = 0, x1 = 0, y0 = 0, y1 = 0;
        if (pt0.x < pt1.x) {
            x0 = pt0.x;
            y0 = pt0.y;
            x1 = pt1.x;
            y1 = pt1.y;
        } else {
            x0 = pt1.x;
            y0 = pt1.y;
            x1 = pt0.x;
            y1 = pt0.y;
        }
        int dx = x1 - x0;
        int dy = y1 - y0;
        int x = 0, y = 0;
        int maxi = 0, mini = 0, maxc = 0;
        maxc = (abs(dx) > abs(dy)) ? abs(dx) : abs(dy);
        int ind = 0;//ÒòÎªÏÂ±ê´Ó0¿ªÊ¼£»
        uchar *dataPtr;
        Mat lineMask(maxc + 1, 2, CV_8U, Scalar(0));
        if (dx == 0) {
            x = x0;
            //if (dy < 0)
            {
                maxi = (y0 > y1) ? y0 : y1;
                mini = (y0 < y1) ? y0 : y1;
            }
            for (y = mini; y <= maxi; y++) {
                dataPtr = lineMask.ptr<uchar>(ind);
                dataPtr[0] = (uchar) y;
                dataPtr[1] = (uchar) x;
                ind++;
            }
        } else if (abs(dy) > abs(dx)) {
            //if (dy < 0)
            {
                maxi = (y0 > y1) ? y0 : y1;
                mini = (y0 < y1) ? y0 : y1;
            }
            for (y = mini; y <= maxi; y++) {
                x = (int) (((dx + 0.0000) / dy + 0.0000) * (y - y0) + x0);
                dataPtr = lineMask.ptr<uchar>(ind);
                dataPtr[0] = (uchar) y;
                dataPtr[1] = (uchar) x;
                ind++;
            }
        } else {
            for (x = x0; x <= x1; x++) {
                y = (int) (((dy + 0.0000) / dx + 0.0000) * (x - x0) + y0);
                dataPtr = lineMask.ptr<uchar>(ind);
                dataPtr[0] = (uchar) y;
                dataPtr[1] = (uchar) x;
                ind++;
            }
        }
        return lineMask;
    }

    Mat drawLine(Point pt0, Point pt1, Mat orgmask) {
        Mat mask = orgmask.clone();
        int maxc = 0;
        int dx = pt0.x - pt1.x;
        int dy = pt0.y - pt1.y;
        maxc = (abs(dx) > abs(dy)) ? abs(dx) : abs(dy);
        Mat tmpMat(maxc + 1, 2, CV_8U, Scalar(0));
        int tmp = 0, tmp1 = 0;
        tmpMat = getLinePts(pt0, pt1);
        uchar *dataPtr;
        uchar *dataPtr1;
        for (int i = 0; i < maxc + 1; i++) {
            dataPtr1 = tmpMat.ptr<uchar>(i);
            tmp = dataPtr1[0];
            dataPtr = mask.ptr<uchar>(tmp);
            tmp1 = dataPtr1[1];
            dataPtr[tmp1] = 255;
        }

        return mask;
    }


    Mat getBaseMask(int theta, int masksize) {
        Mat mask(masksize, masksize, CV_8U, Scalar(0));
        int halfsize = (masksize - 1) / 2;
        int x0 = 0, y0 = 0, x1 = 0, y1 = 0;
        Point pt0, pt1;
        uchar *dataPtr;
        if (theta == 0) {
            dataPtr = mask.ptr<uchar>(halfsize);
            for (int i = 0; i < mask.cols; i++)
                dataPtr[i] = 255;
        } else if (theta == 90) {
            for (int i = 0; i < mask.rows; i++) {
                dataPtr = mask.ptr<uchar>(i);
                dataPtr[halfsize] = 255;
            }
        } else {
            x0 = -halfsize;
            y0 = (int) (x0 * (tan(((float) theta + 0.0000) / 180 * PI) + 0.0000));
            if (y0 < -halfsize) {
                y0 = -halfsize;
                x0 = (int) (y0 * ((1 / (tan((float) theta * PI / 180) + 0.0000) + 0.0000)));
            }
            x1 = halfsize;
            y1 = (int) (x1 * (tan((float) theta * PI / 180) + 0.0000));
            if (y1 > halfsize) {
                y1 = halfsize;
                x1 = (int) (y1 * ((1 / (tan((float) theta * PI / 180) + 0.0000) + 0.0000)));
            }

            pt0.y = halfsize - y0;
            pt0.x = halfsize + x0;
            pt1.y = halfsize - y1;
            pt1.x = halfsize + x1;
            mask = drawLine(pt0, pt1, mask);
        }
        return mask;
    }

    Mat rotatex(Mat mask) {
        int h = 0, w = 0;
        h = mask.rows;
        w = mask.cols;
        Mat rotatedMask(h, w, CV_8U, Scalar(0));
        uchar *dataPtr, *dataPtr2;
        for (int i = 0; i < h; i++) {
            for (int j = 0; j < w; j++) {
                dataPtr = mask.ptr<uchar>(i);
                dataPtr2 = rotatedMask.ptr<uchar>(i);
                dataPtr2[j] = dataPtr[w - j - 1];
            }
        }
        return rotatedMask;
    }


    Mat getLineMask(int theta, int masksize) {
        Mat mask;
        Mat linemask;
        if (theta > 90) {
            mask = getBaseMask(180 - theta, masksize);
            linemask = rotatex(mask);
        } else {
            linemask = getBaseMask(theta, masksize);
        }
        return linemask;
    }

    Mat myIm2Double(Mat in) {
        Mat result;
        if (in.depth() == CV_32F) {
            return in;
        }
        if (in.depth() == CV_32F) {
            in.convertTo(result, CV_32FC(in.channels()), 1.0);
            return result;
        }
        CHECK_EQ(in.depth(), CV_8U);
        in.convertTo(result, CV_32FC(in.channels()), 1.0 / 255.0);
        return result;
    }

/*
 * Imput image should be a 3-channel image
 */
    Mat create_mask(Mat _img, double thresh) {
        CHECK(_img.channels() == 3);
        Mat img = myIm2Double(_img);
        vector<Mat> image_channels;
        split(img, image_channels);
        Mat rChannel(_img.rows, _img.cols, CV_32FC1, Scalar(0));
        for (int i = 0; i < 3; i++) {
            rChannel += image_channels[i].mul(image_channels[i]);
        }
        double maxValue, minValue;
        minMaxLoc(rChannel, &minValue, &maxValue);
        Mat seed = rChannel < thresh * maxValue;
        seed = 255 - seed;
        Mat blur;
        medianBlur(seed, blur, 11);
        return blur;
    }

////////////////////////////////////////////////////////////////////////////
///////////////////////////////////Peak finder//////////////////////////////
    Mat findPeak(Mat, Mat);

#include <limits.h>
#include <math.h>
#include <stdbool.h>
#include <string.h>

//VXparam_t par[] = {  // prefix, value,  description
//    {"if=", 0, " input image"},
//    {0, 0, 0}};

//#define IVAL par[0].val

#define ARRAY_WIDTH 3
#define MAX_LMR 60000
#define ARRAYSIZE(x) (sizeof(x) / sizeof(*(x)))

#define MIN_DIFF 2
//#define MIN_HEIGHT 3
#define MAX_GAP 3
#define THRESH_TOP 1  // the intensity range of the top
#define SLOPE_MAX 500
#define ELEMENT 31  // number of elements in array p[]
#define NUM_PRO 30  // corresponds to 6 degree cross profile
//#define I 11        // for test

    /**
     *
     * @param im inverse green level image
     * @param _lmr lmr information, arrays of row_num, col_num and label
     * @return a list of peaks and their related information
     */
    Mat findPeak(Mat im, Mat _lmr) {
        //  VXparse(&argc, &argv, par);
        //  Vfstruct(im);
        //  Vfread(&im, IVAL);s
        assert(sizeof(int) == 4);
        assert(im.channels() == 1);
        assert(_lmr.depth() == CV_32S);
//        assert(false);

        Mat array;
        _lmr.convertTo(array, CV_32S);

        int p[NUM_PRO][ELEMENT];

        int count_pro = 0, count_ele = 0;

        int max[NUM_PRO], max_d[NUM_PRO];
        int inc_s[NUM_PRO], inc_e[NUM_PRO], dec_s[NUM_PRO], dec_e[NUM_PRO],
                center[NUM_PRO];  // pixel indexes
        int gap[NUM_PRO][ELEMENT], inc[NUM_PRO][ELEMENT], count_gap[NUM_PRO],
                count_inc[NUM_PRO], stop[NUM_PRO];
        int gap_dec[NUM_PRO][ELEMENT], dec[NUM_PRO][ELEMENT], count_gapd[NUM_PRO],
                count_dec[NUM_PRO];

        int xoffset[NUM_PRO][ELEMENT];
        int yoffset[NUM_PRO][ELEMENT];

        long array_in = 0;
        long label_value = 0;
        int y_in = 0, x_in = 0;
        int xx = 0, yy = 0;  // count element in angle;
        double theta = 0;
        int k = 0, i = 0, j = 0;

        // VXparse(&argc,&argv,par);
        long count_out = 0;
        struct Outdata {
            int x_out;
            int y_out;
            long label;
            float mpwidths;
            float sdpwidths;
            float mtwidths;
            float sdtwidths;
            float sdrslopes;
            float cvrheights;
            float cvpheights;
            float score;
        };
        struct Outdata Out[MAX_LMR];

        //  long array[MAX_LMR][ARRAY_WIDTH];
        //
        //  FILE *fp_in;
        //
        //  fp_in = fopen("lmrout.csv", "r");
        //  char buffer[MAX_LMR], *ptr;
        //  size_t read_i, read_j, read_k;
        //
        //  if (!fp_in) {
        //    fprintf(stderr, "failed to open local maximum record file!\n");
        //  } else {
        //    for (read_i = 0; fgets(buffer, sizeof buffer, fp_in); read_i++) {
        //      for (read_j = 0, ptr = buffer; read_j < ARRAYSIZE(*array);
        //           ++read_j, ++ptr) {
        //        array[read_i][read_j] = strtol(ptr, &ptr, 10);
        //        fprintf(stderr, "array[%ld][%ld]=%ld\n", read_i, read_j,
        //                array[read_i][read_j]);
        //      }
        //    }
        //  }
        //
        //  fclose(fp_in);

        size_t read_i = (size_t) array.rows;

//        fprintf(stderr, "read_i=%ld\n", read_i);

        // start read input csv file
        // input y,x axises, and corresponding label value
        for (array_in = 0; array_in < read_i; array_in++) {
//            fprintf(stderr, "\ncurrently processing No.%ld cell\n", array_in);
            y_in = (int) array.at<int>(array_in, 0);
            x_in = (int) array.at<int>(array_in, 1);
            label_value = array.at<int>(array_in, 2);
//            fprintf(stderr, "y,x,label=%d %d %d\n", y_in, x_in, label_value);
            if (y_in >= im.rows - 15 || y_in <= 0 + 15 || x_in >= im.cols - 15 ||
                x_in <= 0 + 15)
                continue;

            // y_in=atoi(Y);
            // x_in=atoi(X);

            // fprintf(stderr,"output profile\n");
            // profile scanning
            // TODO Move this to outer loop
            for (k = 0; k < NUM_PRO; k++) {
                theta = 6 * k;  // angle
                if (theta < 45 || theta > 135) {
                    i = 0;
                    for (xx = -15; xx <= 15; xx++) {
                        yy = (int) round(xx * tan(theta / 180.0 * 3.1415926));
                        xoffset[k][i] = xx;
                        yoffset[k][i] = yy;
                        i++;
                    }
                } else {
                    i = 0;
                    for (yy = -15; yy <= 15; yy++) {
                        xx = (int) round(yy / tan(theta / 180.0 * 3.1415926));
                        xoffset[k][i] = xx;
                        yoffset[k][i] = yy;
                        i++;
                    }
                }
            }

            // Get their pixel values
            for (k = 0; k < NUM_PRO; k++) {
                int xx, yy;
                for (i = 0; i < ELEMENT; i++) {
                    xx = xoffset[k][i];
                    yy = yoffset[k][i];
                    p[k][i] = im.at<uchar>(y_in + yy, x_in + xx);
//                    fprintf(stderr , "val %d \n",(int) im.at<uchar>(y_in + yy, x_in + xx));
                }
            }

            for (i = 0; i < NUM_PRO; i++)
                for (j = 0; j < ELEMENT; j++) {
                    gap[i][j] = 0;
                    inc[i][j] = 0;
                    gap_dec[i][j] = 0;
                    dec[i][j] = 0;
                }

            for (i = 0; i < NUM_PRO; i++) {
                max[i] = 0;
                max_d[i] = 0;
                inc_s[i] = 0;
                inc_e[i] = ELEMENT / 2;
                dec_s[i] = ELEMENT / 2 + 1;
                dec_e[i] = ELEMENT - 1;
                center[i] = ELEMENT / 2;
                stop[i] = 0;
                count_gap[i] = 0;
                count_inc[i] = 0;
                count_gapd[i] = 0;
                count_dec[i] = 0;
            }

            for (count_pro = 0; count_pro < NUM_PRO; count_pro++) {
                for (count_ele = 0; count_ele < ELEMENT; count_ele++) {  // peak detection
                    if (p[count_pro][count_ele] > max[count_pro]) {
                        max[count_pro] = p[count_pro][count_ele];  // record peak value
                        inc_e[count_pro] = count_ele;  // record peak index, also the end
                        // index of rising ramp
                    }
                }
            }

            // find starting point of decreasing ramp
            for (count_pro = 0; count_pro < NUM_PRO; count_pro++) {
                for (count_ele = inc_e[count_pro]; count_ele < ELEMENT; count_ele++) {
                    if (max[count_pro] - p[count_pro][count_ele] > THRESH_TOP) {
                        max_d[count_pro] = p[count_pro][count_ele - 1];
                        dec_s[count_pro] = count_ele - 1;
                        break;
                    }
                }
            }

            // calculate elements in gap[]   (increasing)
            for (count_pro = 0; count_pro < NUM_PRO; count_pro++) {
                for (j = inc_e[count_pro]; j > 0; j--) {
                    if (p[count_pro][j] - p[count_pro][j - 1] >= MIN_DIFF) {
                        gap[count_pro][count_gap[count_pro]] = j;
                        count_gap[count_pro]++;
                    } else {
                        if (p[count_pro][j + 1] - p[count_pro][j] >= MIN_DIFF) {
                            gap[count_pro][count_gap[count_pro]] = j;
                            count_gap[count_pro]++;
                        }
                    }
                }
            }

            if (p[count_pro][1] - p[count_pro][0] >= MIN_DIFF) {
                gap[count_pro][count_gap[count_pro]++] = 0;
            }

            // calculate elements in gap_dec[]  (decreasing)
            for (count_pro = 0; count_pro < NUM_PRO; count_pro++) {
                for (j = dec_s[count_pro]; j < ELEMENT - 1; j++) {
                    if (p[count_pro][j] - p[count_pro][j + 1] >= MIN_DIFF) {
                        gap_dec[count_pro][count_gapd[count_pro]] = j;
                        count_gapd[count_pro]++;
                    } else {
                        if (p[count_pro][j - 1] - p[count_pro][j] >= MIN_DIFF) {
                            gap_dec[count_pro][count_gapd[count_pro]] = j;
                            count_gapd[count_pro]++;
                        }
                    }
                }
            }

            if (p[count_pro][ELEMENT - 2] - p[count_pro][ELEMENT - 1] >= MIN_DIFF) {
                gap_dec[count_pro][count_gapd[count_pro]++] = ELEMENT - 1;
            }

            // store increasing ramp's indexes in inc[]
            for (count_pro = 0; count_pro < NUM_PRO; count_pro++) {
                for (j = 0; j < count_gap[count_pro]; j++) {
                    if (gap[count_pro][j] - gap[count_pro][j + 1] > MAX_GAP) {
                        if (gap[count_pro][j - 1] - gap[count_pro][j] <= MAX_GAP) {
                            inc[count_pro][j] = gap[count_pro][j];
                            count_inc[count_pro]++;
                        }
                        break;
                    } else {
                        inc[count_pro][j] = gap[count_pro][j];
                        count_inc[count_pro]++;
                    }
                }
            }

            for (count_pro = 0; count_pro < NUM_PRO; count_pro++) {
                inc_s[count_pro] = inc[count_pro][count_inc[count_pro] - 1];
            }

            // store decreasing ramp's indexes in dec[]
            for (count_pro = 0; count_pro < NUM_PRO; count_pro++) {
                for (j = 0; j < count_gapd[count_pro]; j++) {
                    if (abs(gap_dec[count_pro][j] - gap_dec[count_pro][j + 1]) > MAX_GAP) {
                        if (abs(gap_dec[count_pro][j - 1] - gap_dec[count_pro][j]) <=
                            MAX_GAP) {
                            dec[count_pro][j] = gap_dec[count_pro][j];
                            count_dec[count_pro]++;
                        }
                        break;
                    } else {
                        dec[count_pro][j] = gap_dec[count_pro][j];
                        count_dec[count_pro]++;
                    }
                }
            }

            if (gap_dec[count_pro][count_gapd[count_pro] - 1] -
                gap_dec[count_pro][count_gapd[count_pro] - 2] >
                MAX_GAP) {
            }

            for (count_pro = 0; count_pro < NUM_PRO; count_pro++) {
                if (dec[count_pro][count_dec[count_pro] - 1] != 0) {
                    dec_e[count_pro] = dec[count_pro][count_dec[count_pro] - 1];
                }
            }

            // calculate profile property
            int w_peak[NUM_PRO], w_top[NUM_PRO], h_inc[NUM_PRO], h_dec[NUM_PRO],
                    h_peak[NUM_PRO];                   // profile properties
            float s_inc[NUM_PRO], s_dec[NUM_PRO];  // increasing and decreasing slope

            for (i = 0; i < NUM_PRO; i++) {
                w_peak[i] = 0;
                w_top[i] = 1;
                h_inc[i] = 0;
                h_dec[i] = 0;
                h_peak[i] = 0;
                s_inc[i] = 0;
                s_dec[i] = 0;
            }

            // if only exists rising ramp, deceasing index starts and ends at last index
            for (count_pro = 0; count_pro < NUM_PRO; count_pro++) {
                if (inc_e[count_pro] == ELEMENT - 1) {
                    dec_s[count_pro] = ELEMENT - 1;
                    dec_e[count_pro] = ELEMENT - 1;
                }
            }
            // if only exists descending ramp, increasing index starts and ends at first
            // index
            for (count_pro = 0; count_pro < NUM_PRO; count_pro++) {
                if (dec_s[count_pro] == 0) {
                    inc_s[count_pro] = 0;
                    inc_e[count_pro] = 0;
                }
            }

            for (count_pro = 0; count_pro < NUM_PRO; count_pro++) {
                if (dec_s[count_pro] != inc_e[count_pro]) {
                    w_top[count_pro] = dec_s[count_pro] - inc_e[count_pro];
                } else {
                    w_top[count_pro] = 1;
                }
                w_peak[count_pro] = dec_e[count_pro] - inc_s[count_pro];
                h_inc[count_pro] =
                        p[count_pro][inc_e[count_pro]] - p[count_pro][inc_s[count_pro]];
                h_dec[count_pro] =
                        p[count_pro][dec_s[count_pro]] - p[count_pro][dec_e[count_pro]];
                center[count_pro] = (inc_e[count_pro] + dec_s[count_pro]) / 2;
            }

            // calculate increasing slope (rise&descend / only descending ramp)
            for (count_pro = 0; count_pro < NUM_PRO; count_pro++) {
                if (inc_s[count_pro] != inc_e[count_pro]) {
                    s_inc[count_pro] =
                            (float) (h_inc[count_pro] / (inc_e[count_pro] - inc_s[count_pro]));
                } else {
                    s_inc[count_pro] = SLOPE_MAX;
                }
            }

            // calculate decreasing slope (rise&descend / only rising ramp)
            for (count_pro = 0; count_pro < NUM_PRO; count_pro++) {
                if (dec_s[count_pro] != dec_e[count_pro]) {
                    s_dec[count_pro] =
                            (float) (h_dec[count_pro] / (dec_e[count_pro] - dec_s[count_pro]));
                } else {
                    s_dec[count_pro] =
                            SLOPE_MAX;  // decreasing slope is also considered as positive
                }
            }

            // calculate peak height
            for (count_pro = 0; count_pro < NUM_PRO; count_pro++) {
                if (p[count_pro][inc_s[count_pro]] > p[count_pro][dec_e[count_pro]]) {
                    h_peak[count_pro] =
                            p[count_pro][center[count_pro]] - p[count_pro][inc_s[count_pro]] +
                            ((p[count_pro][inc_s[count_pro]] - p[count_pro][dec_e[count_pro]]) *
                             (center[count_pro] - inc_s[count_pro]) / w_peak[count_pro]);
                } else if (p[inc_s[count_pro]] < p[dec_e[count_pro]]) {
                    h_peak[count_pro] =
                            p[count_pro][center[count_pro]] - p[count_pro][dec_e[count_pro]] +
                            ((p[count_pro][dec_e[count_pro]] - p[count_pro][inc_s[count_pro]]) *
                             (dec_e[count_pro] - center[count_pro]) / w_peak[count_pro]);
                } else {
                    h_peak[count_pro] =
                            p[count_pro][center[count_pro]] - p[count_pro][inc_s[count_pro]];
                }
            }

            for (count_pro = 0; count_pro < NUM_PRO; count_pro++) {
                if (!w_top[count_pro] && !w_peak[count_pro] && !h_inc[count_pro] &&
                    !h_dec[count_pro]) {
                    s_inc[count_pro] = 0;
                    s_dec[count_pro] = ELEMENT - 1;
                    h_peak[count_pro] = 0;
                    center[count_pro] = ELEMENT / 2;
                }
            }

            for (count_pro = 0; count_pro < NUM_PRO; count_pro++) {
            }
//            fprintf(stderr, "start calculating properties...\n");

            float twidths[NUM_PRO];       // top width
            float pwidths[NUM_PRO];       // peak width
            float rheights[NUM_PRO * 2];  // ramp height
            float pheights[NUM_PRO];      // peak height
            float rslopes[NUM_PRO * 2];   // ramp slope

            int ct = 0;
            for (count_pro = 0; count_pro < NUM_PRO; count_pro++) {
                if (count_inc[count_pro] && count_dec[count_pro]) {
                    twidths[ct] = w_top[count_pro];
                    pwidths[ct] = w_peak[count_pro];
                    rheights[2 * ct] = h_inc[count_pro];
                    rheights[2 * ct + 1] = h_dec[count_pro];
                    rslopes[2 * ct] = s_inc[count_pro];
                    rslopes[2 * ct + 1] = s_dec[count_pro];
                    pheights[ct] = h_peak[count_pro];
                    ct++;
                }
            }

            // calculate profile property statistics
            float mpwidths, mtwidths, mrslopes, mrheights, mpheights, MINpheights;
            float sdpwidths, sdtwidths, sdrslopes, sdrheights, sdpheights;
            float cvrheights, cvpheights;
            float score;

            // twidths
            float sum1 = 0, var1 = 0;  // calculate mean twidths
            for (i = 0; i < ct; i++) {
                sum1 += twidths[i];
            }
            mtwidths = sum1 / (float) ct;

            for (j = 0; j < ct; j++) {
                var1 += (twidths[j] - mtwidths) * (twidths[j] - mtwidths);
            }
            sdtwidths = sqrt(var1 / (float) ct);

            // pwidths
            float sum2 = 0, var2 = 0;  // calculate mean pwidths
            for (i = 0; i < ct; i++) {
                sum2 += pwidths[i];
            }
            mpwidths = sum2 / (float) ct;

            for (j = 0; j < ct; j++) {
                var2 += (pwidths[j] - mpwidths) * (pwidths[j] - mpwidths);
            }
            sdpwidths = sqrt(var2 / (float) ct);

            // rslopes
            float sum3 = 0, var3 = 0;  // calculate mean rslopes
            for (i = 0; i < 2 * ct; i++) {
                sum3 += rslopes[i];
            }
            mrslopes = sum3 / (2 * (float) ct);

            for (j = 0; j < 2 * ct; j++) {
                var3 += (rslopes[j] - mrslopes) * (rslopes[j] - mrslopes);
            }
            sdrslopes = sqrt(var3 / (2 * (float) ct));

            // rheights
            float sum4 = 0, var4 = 0;  // calculate mean mrheights
            for (i = 0; i < 2 * ct; i++) {
                sum4 += rheights[i];
            }
            mrheights = sum4 / (2 * (float) ct);

            for (j = 0; j < 2 * ct; j++) {
                var4 += (rheights[j] - mrheights) * (rheights[j] - mrheights);
            }
            sdrheights = sqrt(var4 / (2 * (float) ct));

            cvrheights = sdrheights / mrheights;

            // pheights
            float sum5 = 0, var5 = 0;  // calculate mean mpheights
            for (i = 0; i < ct; i++) {
                sum5 += pheights[i];
            }
            mpheights = sum5 / (float) ct;

            for (j = 0; j < ct; j++) {
                var5 += (pheights[j] - mpheights) * (pheights[j] - mpheights);
            }
            sdpheights = sqrt(var5 / (float) ct);

            cvpheights = sdpheights / mpheights;

            MINpheights = pheights[0];
            for (i = 0; i < ct; i++) {
                if (pheights[i] < MINpheights) {
                    MINpheights = pheights[i];
                }
            }

            // score
            score = (MINpheights * mrslopes) /
                    (1 + sdpwidths + sdtwidths + sdrslopes + sdrheights + sdpheights);

            if (mpwidths > INT_MIN && mtwidths > INT_MIN && mpheights > INT_MIN &&
                cvpheights > INT_MIN && cvpheights < INT_MAX && score >= 0) {
                Out[count_out].y_out = y_in;
                Out[count_out].x_out = x_in;
                Out[count_out].mpwidths = mpwidths;
                Out[count_out].sdpwidths = sdpwidths;
                Out[count_out].mtwidths = mtwidths;
                Out[count_out].sdtwidths = sdtwidths;
                Out[count_out].sdrslopes = sdrslopes;
                Out[count_out].cvrheights = cvrheights;
                Out[count_out].cvpheights = cvpheights;
                Out[count_out].score = score;
                Out[count_out].label = label_value;

                count_out++;

                // fprintf(stderr, "%d,%d,%ld,%f,%f,%f,%f,%f,%f,%f,%f\n", y_in, x_in,
                //         label_value, mpwidths, sdpwidths, mtwidths, sdtwidths, sdrslopes,
                //         cvrheights, cvpheights, score);
            }
        }

//        FILE *fp_out;
//        fp_out = fopen("peakout.csv", "w");
        Mat result = cv::Mat((int) count_out, 11, CV_32F);

        for (i = 0; i < count_out; i++) {
//            fprintf(stderr, "write %d,%d,%ld,%f,%f,%f,%f,%f,%f,%f,%f\n", Out[i].y_out,
//                    Out[i].x_out, Out[i].label, Out[i].mpwidths, Out[i].sdpwidths,
//                    Out[i].mtwidths, Out[i].sdtwidths, Out[i].sdrslopes,
//                    Out[i].cvrheights, Out[i].cvpheights, Out[i].score);
//            fprintf(fp_out, "%d,%d,%ld,%f,%f,%f,%f,%f,%f,%f,%f\n", Out[i].y_out,
//                    Out[i].x_out, Out[i].label, Out[i].mpwidths, Out[i].sdpwidths,
//                    Out[i].mtwidths, Out[i].sdtwidths, Out[i].sdrslopes,
//                    Out[i].cvrheights, Out[i].cvpheights, Out[i].score);
            result.at<float>(i, 0) = Out[i].y_out;
            result.at<float>(i, 1) = Out[i].x_out;
            result.at<float>(i, 2) = Out[i].label;
            result.at<float>(i, 3) = Out[i].mpwidths;
            result.at<float>(i, 4) = Out[i].sdpwidths;
            result.at<float>(i, 5) = Out[i].mtwidths;
            result.at<float>(i, 6) = Out[i].sdtwidths;
            result.at<float>(i, 7) = Out[i].sdrslopes;
            result.at<float>(i, 8) = Out[i].cvrheights;
            result.at<float>(i, 9) = Out[i].cvpheights;
            result.at<float>(i, 10) = Out[i].score;
        }

//        fclose(fp_out);

        // exit(0);
        return result;
    }

////////////////////////////////////////////////////////////////////////////
    Mat getMask(Mat _img, double thresh = 0.0001) {
        Mat img;
        _img.convertTo(img, CV_32FC(_img.channels()), 1.0 / 255.0);
        vector<Mat> image_channels;
        split(img, image_channels);
        Mat rChannel(img.rows, img.cols, CV_32FC1, Scalar(0));
        for (int i = 1; i < 2; i++) {
            rChannel += image_channels[i].mul(image_channels[i]);
        }
        double maxValue, minValue;
        minMaxLoc(rChannel, &minValue, &maxValue);
        return rChannel > thresh * maxValue;
    }

    Mat exudate_extract(Mat image) {
        //Mat image = imread("im01.jpg");//¶ÁÈëÔ­Ê¼Í¼Ïñ
        //Mat mask = imread("mask.tif",0);//ÑÚÄ¤
        Mat mask = getMask(image);
        Mat image_small = mask.clone();
        resize(image, image_small, image_small.size());//µ÷ÕûÍ¼Ïñ´óÐ¡ÒÔÊÊÓ¦ÑÚÄ¤
        image = image_small.clone();


        //½«Ô­Ê¼Í¼Ïñ×ª»»Îª»Ò¶ÈÍ¼ÏñÒÔ±¸ºóÐø´¦Àí
        Mat grey;
        cvtColor(image, grey, CV_BGR2GRAY);

        //±Õ²Ù×÷³ýÈ¥Ñª¹Ü£¬±£´æÔÚe1ÖÐ
        Mat e1;
        Mat element1 = getStructuringElement(MORPH_ELLIPSE, Size(7, 7));
        morphologyEx(grey, e1, CV_MOP_CLOSE, element1);
//        LOG(INFO) << "1";
        /*ÇóÈ¡Éø³ö±ß½ç*/
        int m = e1.rows, n = e1.cols;
        Mat e2 = Mat::zeros(m, n, CV_8UC1);//e2µÄ³õÊ¼»¯£¬ÎªÈ«ºÚÄ£°å
        Mat window = Mat::zeros(11, 11, CV_8UC1);//¶¨ÒåÒ»¸ö11*11µÄ¾í»ýµ¥Ôª
        int flag;
        for (int i = 5; i <= m - 6; ++i)
            for (int j = 5; j <= n - 6; ++j) {
                flag = true;
                for (int a = 0; a < 11; ++a)
                    for (int b = 0; b < 11; ++b) {
                        window.at<uchar>(a, b) = e1.at<uchar>(i - 5 + a,
                                                              j - 5 + b);//window±£´æ´Óe1ÖÐÈ¡³öµÄ11*11µÄÐ¡Í¼¿é
                        if (mask.at<uchar>(i - 5 + a, j - 5 + b) != 0)flag = false;
                    }
                if (flag)continue;

                Scalar mean, stdDev;
                meanStdDev(window, mean, stdDev, noArray());//¼ÆËã·½²î£¬²¢±£´æÔÚe2ÖÐ
                e2.at<uchar>(i, j) = stdDev.val[0];
            }

//        LOG(INFO) << "2";

        /*¶Ôe2ãÐÖµ·Ö¸î£¬µÃµ½Éø³öÇøÓòµÄ±ß½ç£¬±£´æÔÚe3ÖÐ*/
        Mat e3 = e2.clone();
        for (int i = 0; i < m; ++i)
            for (int j = 0; j < n; ++j) {
                if (mask.at<uchar>(i, j) == 0)
                    e3.at<uchar>(i, j) = 0;
                else if (e2.at<uchar>(i, j) > 5)//5ÎªãÐÖµ£¬¾­ÑéÖµ£¬¿É¾ÝÍ¼µ÷Õû
                    e3.at<uchar>(i, j) = 255;
                else e3.at<uchar>(i, j) = 0;
            }

        /*¶Ôe3Ö´ÐÐ±Õ²Ù×÷£¬±£Ö¤Éø³öÇøÓòµÄÁ¬Í¨ÐÔ£¬±£´æÔÚe4ÖÐ*/
        Mat e4;
        Mat element2 = getStructuringElement(MORPH_ELLIPSE, Size(13, 13));
        morphologyEx(e3, e4, CV_MOP_CLOSE, element2);

        /*ÔÚÔ­Í¼£¨ÓÃe6¸´ÖÆ£©ÖÐ±ê³öÉø³öÇøÓò*/
        Mat e6 = e1.clone();
        for (int i = 0; i < m; ++i)
            for (int j = 0; j < n; ++j) {
                if (e4.at<uchar>(i, j) != 0)
                    e6.at<uchar>(i, j) = 0;
            }

//        LOG(INFO) << "3";

        Mat element3 = getStructuringElement(MORPH_ELLIPSE, Size(2, 2));//ÐÎÌ¬Ñ§²Ù×÷µÄºË
        Mat e7 = e6.clone();
//        imwrite("out/e6.png", e6);
        Mat e7p;
        int pp = 1;//Ñ­»·±êÖ¾
        /*µü´úÅòÕÍµÃµ½Î´²¡±äÍ¼Ïñ*/
        int prev_change = 999999999;
        int this_change = 999999998;
        while (prev_change - this_change != 0) {
            dilate(e7, e7p, element3);
            prev_change = this_change;
            this_change = 0;
            for (int i = 0; i < m; ++i)
                for (int j = 0; j < n; ++j) {
                    if (e7.at<uchar>(i, j) == 0) {
                        e7.at<uchar>(i, j) = e7p.at<uchar>(i, j);
                        this_change += 1;
                    }
                }
//            LOG(INFO) << "3.5 " << this_change;
        }
//        imwrite("out/e7.png", e7);
//        LOG(INFO) << "4";

        Mat e8 = e1 - e7;//Ô­Ê¼Í¼Ïñ¼õÈ¥Îª²¡±äÍ¼Ïñ¼´Îª²¡±äÇøÓò£¨±£´æÔÚe8ÖÐ£©

        double max_e8;//×îÁÁµã
        minMaxLoc(e8, NULL, &max_e8, NULL, NULL);
        double alpha = 0.15;//alpha=0.15Îª¾­ÑéÖµ£¬¿É¸ù¾ÝÍ¼ÏñÁÁ¶Èµ÷Õû
        int area = 0;
        for (int i = 0; i < m; ++i)
            for (int j = 0; j < n; ++j) {
                if (e8.at<uchar>(i, j) > alpha * max_e8) {
                    image.at<Vec3b>(i, j)[0] = 0;
                    image.at<Vec3b>(i, j)[1] = 0;
                    image.at<Vec3b>(i, j)[2] = 0;
                    ++area;//¼ÇÂ¼Éø³öÃæ»ýµÄÏñËØµã
                }
            }
//        LOG(INFO) << "5";

//        cout << "area = " << area * 1.0 / 146861 << endl;//·ÖÄ¸ÎªÓÐÐ§ÇøÓòÏñËØµãµÄ×ÜÊý£¨¸ù¾ÝmaskµÄ´óÐ¡ºÍÐÎ×´µ÷Õû£©
        return image;
    }

    /////////////////////////////////////////////////////////////////////
    void brightCorrection(Mat &in, Mat &out)
    {
        cvtColor(in, out, CV_BGR2HSV);//½«Ô­À´µÄBGRÍ¨µÀµÄÍ¼Ïñ×ª±äÎªHSVÍ¨µÀµÄÍ¼Ïñ
        for(int y = 0; y < in.rows; y++ ) {
            for(int x = 0; x < in.cols; x++ ) {
                int bright = out.at<Vec3b>(y,x)[2];
                out.at<Vec3b>(y,x)[2]= saturate_cast<uchar>(pow((65025-(255-bright)*(255-bright)),0.5)); //¶ÔVÍ¨µÀ½øÐÐ±ä»»
            }
        }
        cvtColor(out, out, CV_HSV2BGR);//½«HSVÍ¨µÀ»¹Ô­ÎªBGRÍ¨µÀ
    }


/*
¹¦ÄÜ£º¶Ô±È¶ÈÔöÇ¿
Ê¹ÓÃ£ºimageAdjust£¨Ô­Ê¼Í¼Ïñ£¬×îÖÕÍ¼Ïñ£¬ÊäÈëµÍãÐÖµ£¬ÊäÈë¸ßãÐÖµ£¬Êä³öµÍãÐÖµ£¬ Êä³ö¸ßãÐÖµ£¬Ö¸Êý²ÎÊý£©
*/
    int imageAdjust(Mat &src,Mat &dst,  double low_in,double high_in,  double low_out,double high_out,  double gamma)
    {
        if ( low_in < 0 && low_in >1 && high_in <0 && high_in >1
             && low_out < 0 && low_out > 1 && high_out < 0 && high_out > 1
             && low_out > high_out )
            return -1;

        double low2 = low_in*255;
        double high2 = high_in*255;
        double bottom2 = low_out*255;
        double top2 = high_out*255;
        double err_in = high2 - low2;
        double err_out = top2 - bottom2;

        int x, y;
        double val;
        uchar* p;

        for (y = 0; y < src.rows; y++) {
            p = src.ptr<uchar>(y);  //ÐÐÖ¸Õë
            for (x = 0; x < src.cols; x++) {
                val = pow((p[x] - low2)/err_in,gamma)*err_out + bottom2;
                if (val > 255 )  val = 255;
                if (val < 0 )    val = 0;
                dst.at<uchar>(y,x) = saturate_cast<uchar>(val);
            }
        }

    }


/*
¹¦ÄÜ£ºÌáÈ¡Ñª¹ÜÇøÓò
Ê¹ÓÃ£ºvesselExtract£¨ÊäÈëÍ¼Ïñ£¬Êä³öÍ¼Ïñ£¬ãÐÖµ£©
*/
    void vesselExtract(Mat &in, Mat &out, int threshold)
    {
        //KirschËã×Ó£¨h1~h8£©
        //t1~t8ÎªÖ´ÐÐËã×ÓºóµÄ½á¹ûÍ¼Ïñ
        Mat h1=(Mat_<double>(3,3)<<5,-3,-3,5,0,-3, 5,-3,-3);
        h1 /=15;
        Mat t1;
        filter2D (in,t1,in.depth (),h1);

        Mat h2=(Mat_<double>(3,3)<<-3,-3,5,-3,0,5,-3,-3,5);
        h2 /=15;
        Mat t2;
        filter2D (in,t2,in.depth (),h2);

        Mat h3=(Mat_<double>(3,3)<<-3,-3,-3,5,0,-3, 5,5,-3);
        h3 /=15;
        Mat t3;
        filter2D (in,t3,in.depth (),h3);

        Mat h4=(Mat_<double>(3,3)<<-3,5,5,-3,0,5, -3,-3,-3);
        h4 /=15;
        Mat t4;
        filter2D (in,t4,in.depth (),h4);

        Mat h5=(Mat_<double>(3,3)<<-3,-3,-3,-3,0,-3,5,5,5);
        h5 /=15;
        Mat t5;
        filter2D (in,t5,in.depth (),h5);

        Mat h6=(Mat_<double>(3,3)<<5,5,5,-3,0,-3, -3,-3,-3);
        h6 /=15;
        Mat t6;
        filter2D (in,t6,in.depth (),h6);

        Mat h7=(Mat_<double>(3,3)<<-3,-3,-3,-3,0,5, -3,5,5);
        h7 /=15;
        Mat t7;
        filter2D (in,t7,in.depth (),h7);

        Mat h8=(Mat_<double>(3,3)<<5,5,-3,5,0,-3, -3,-3,-3);
        h8 /=15;
        Mat t8;
        filter2D (in,t8,in.depth (),h8);


        out = Mat::zeros(in.rows, in.cols,CV_8UC1);//È«ºÚÄ£°å
        int temp[9] = {0,0,0,0,0,0,0,0,0};

        for(int i=0; i<in.rows; i++)
            for(int j=0; j<in.cols; j++) {
                temp[1] = saturate_cast<uchar>(t1.at<uchar>(i,j));
                temp[2] = saturate_cast<uchar>(t2.at<uchar>(i,j));
                temp[3] = saturate_cast<uchar>(t3.at<uchar>(i, j));
                temp[4] = saturate_cast<uchar>(t4.at<uchar>(i, j));
                temp[5] = saturate_cast<uchar>(t5.at<uchar>(i, j));
                temp[6] = saturate_cast<uchar>(t6.at<uchar>(i, j));
                temp[7] = saturate_cast<uchar>(t7.at<uchar>(i, j));
                temp[8] = saturate_cast<uchar>(t8.at<uchar>(i, j));

                int max=0;
                for(int k=1; k <= 8; ++k)
                    max = (max > temp[k]) ? max : temp[k];//È¡8¸ö¾ØÕóÔËËãºóµÄ×î´óÖµ
                if(max > threshold)out.at<uchar>(i, j) = cv::saturate_cast<uchar>(max);//Èç¹û×î´óÖµ´óÓÚãÐÖµÔò½«¸Ã×î´óÖµ×÷ÎªÊä³ö
            }
    }


Mat hemorrhage_extract(Mat inImg) {
        Mat grey;
        cvtColor(inImg, grey, CV_BGR2GRAY);//½«Ô­Ê¼²ÊÉ«Í¼Ïñ(inImg)×ª»»Îª»Ò¶ÈÍ¼Ïñ£¨grey£©

        /*´´½¨Ô­Ê¼Í¼Ïñ¸±±¾£¬ÒÔ¹©ºóÐøÊ¹ÓÃ*/
        Mat I_initial = inImg.clone();
        Mat I1 = inImg.clone();

        /*Èç¹ûÔ­Ê¼Í¼ÏñÎª²ÊÉ«Í¼Ïñ£¬ÔòÖ´ÐÐÒÔÏÂ²Ù×÷£¨ÁÁ¶ÈÐÞÕýºÍ¶Ô±È¶ÈÔöÇ¿£©£»·ñÔò£¬ÔòÌø¹ý´Ë²½*/
        if (inImg.channels() == 3) {
            /*ÁÁ¶ÈÐÞÕý*/
            Mat BC_inImg = inImg.clone();
            brightCorrection(inImg, BC_inImg);//µ÷ÓÃ¡°ÁÁ¶ÈÐÞÕý¡±º¯Êý
            inImg = BC_inImg.clone();

            /*½«ÁÁ¶ÈÐÞÕýºóµÄ²ÊÉ«Í¼ÏñÍ¼Ïñ×ª±äÎª»Ò¶ÈÍ¼Ïñ*/
            Mat BC_grey;
            cvtColor(inImg, BC_grey, CV_BGR2GRAY);
            inImg = BC_grey.clone();

            /*¶Ô±È¶ÈÔöÇ¿*/
            Mat in_adj = BC_grey.clone();
            imageAdjust(BC_grey, in_adj, 0.3, 0.7, 0, 1,
                        1);//µ÷ÓÃ¡°¶Ô±È¶ÈÔöÇ¿¡±º¯Êý£¬½«0.3-0.7µÄ»Ò¶È·¶Î§À­Éìµ½0-1£¨Ïà¶ÔÓÚ×î´ó»Ò¶ÈÖµ255¶øÑÔ£©£¬0.3ºÍ0.7¾ùÎª¾­ÑéÖµ
            inImg = in_adj.clone();//inImgÖÐ±£´æ¼Ó¹¤ºóµÄ»Ò¶ÈÍ¼
            I1 = in_adj.clone();//I1ÖÐÍ¬Ñù±£´æ¼Ó¹¤ºóµÄ»Ò¶ÈÍ¼
        }

        /*Ñª¹ÜÌáÈ¡*/
        Mat bloodVessels = inImg.clone();
        vesselExtract(inImg, bloodVessels, 18);//µ÷ÓÃ¡°Ñª¹ÜÌáÈ¡¡±º¯Êý£¬18ÎªãÐÖµ£¨¾­ÑéÖµ£¬¿É¾ÝÍ¼ÊÊµ±ÐÞÕý£©

        /*½«ÌáÈ¡µÄÑª¹ÜÇøÓò½øÒ»²½¶þÖµ»¯*/
        for (int y = 0; y < bloodVessels.rows; y++) {
            for (int x = 0; x < bloodVessels.cols; x++) {
                if (bloodVessels.at<uchar>(y, x) > 0)
                    bloodVessels.at<uchar>(y, x) = 255;
            }
        }

        /*¶ÔÌáÈ¡µÄÑª¹Ü²¿Î»Ö´ÐÐ¡°±Õ²Ù×÷¡±£¬ÃÖ²¹Ï¸Ð¡¿ÕÏ¶ÒÔÈ·±£Ñª¹ÜÁ¬Í¨ÐÔ*/
        Mat bV_close;
        Mat element1 = getStructuringElement(MORPH_ELLIPSE, Size(3, 3));
        morphologyEx(bloodVessels, bV_close, CV_MOP_CLOSE, element1);
        bloodVessels = bV_close.clone();

        /*¶ÔÔ­»Ò¶ÈÍ¼ÏñãÐÖµ»¯£¬³õ²½È·¶¨³öÑªÇøÓò*/
        Mat I2;
        threshold(I1, I2, 95, 255, THRESH_BINARY);//135ÎªãÐÖµ£¨¾­ÑéÖµ£¬¿É¾ÝÍ¼µ÷Õû£©£¬255Îª»Ò¶ÈÍ¼ÏñÔÊÐíµÄ×î´óÖµ

        /*È·¶¨ÓÐÐ§¼ÆËãÇøÓò£¬¼´ÔÚ¼ÆËã³öÑª±ÈÀýÊ±µÄ×ÜÃæ»ý*/

//        ofstream Savefile("hemorrhage.txt");
//        Savefile << "³öÑªÇøÓòÏñËØµãÎª£º\n";

        int count = 0, count_blood = 0;
        for (int y = 0; y < grey.rows; y++) {
            for (int x = 0; x < grey.cols; x++) {
                if (grey.at<uchar>(y, x) > 10)//10Îª¾­ÑéÖµ£¬Ö÷ÒªÎªÁËÅÅ³ýÑÛµ×Í¼ÏñÖÜÎ§ºÚÉ«ÂÖÀªµÄ¸ÉÈÅ
                {
                    grey.at<uchar>(y, x) = 255;
                    count++;//¼ÇÂ¼ÓÐÐ§¼ÆËãÇøÓòµÄÏñËØµã
                }
            }
        }

        /*ÔÚ³öÑªÇøÓòÖÐ£¬ÌÞ³ýÑª¹ÜÇøÓò*/
        for (int y = 0; y < I2.rows; y++) {
            for (int x = 0; x < I2.cols; x++) {
                if ((I2.at<uchar>(y, x) == 0) && bloodVessels.at<uchar>(y, x) == 0) {
                    I_initial.at<Vec3b>(y, x)[0] = 0;
                    I_initial.at<Vec3b>(y, x)[1] = 0;
                    I_initial.at<Vec3b>(y, x)[2] = 0;
                    if (grey.at<uchar>(y, x) == 255) {
                        ++count_blood;//¼ÇÂ¼³öÑªÇøÓòµÄÏñËØµã
//                        Savefile << "(" << y << "," << x << ")\n";
                    }
                }
            }
        }

        /*¶Ô³öÑªÇøÓòÖ´ÐÐ£¨Î¢Ð¡µÄ£©¿ª²Ù×÷£¬±£Ö¤³öÑªÇøÓòµÄÁ¬Í¨ÐÔ*/
        Mat result_blood;
        Mat element2 = getStructuringElement(MORPH_ELLIPSE, Size(2, 2));
        morphologyEx(I_initial, result_blood, CV_MOP_OPEN, element2);

        double area;
        area = count_blood * 1.0 / count;//¼ÆËã³öÑª±ÈÀý
//        Savefile << "³öÑª¼ì²â³É¹¦£¡\n";
//        Savefile << "³öÑªÇøÓòÃæ»ý = " << area << "\n";
//        Savefile.close();
//        cout << "³öÑª¼ì²â³É¹¦£¡" << endl;
//        cout << "³öÑªÇøÓòÃæ»ý = " << area << endl;
        return result_blood;
    }



    inline int max(int a, int b) {
        if (a > 0) return a;
        else return b;
    }

    inline int min(int a, int b) {
        if (a < 0) return a;
        else return b;
    }

    Point pol2cart(double a, int b) {
        int x, y;
        x = cvRound(b * cos(a * 0.1 * PI)) + 80;
        y = cvRound(b * sin(a * 0.1 * PI)) + 80;
        return Point(x, y);
    }

    Mat getOpticDisc(Mat img) {
        if (img.empty()) {
//            cout << "Can not open the Image!" << endl << endl;
            return Mat();
        }

        //?????????Green?????

        vector<Mat> channels;
        split(img, channels);
        Mat greenChannel, redChannel;
        greenChannel = channels[1];
        redChannel = channels[2];

        //equalizeHist(greenChannel,greenChannel);imshow("equalization",greenChannel);
        //greenChannel=greenChannel*1.1;

        //??????;
        Mat filtedImage;
        blur(greenChannel, filtedImage, Size(31, 31));
        //imshow("filtedImage",filtedImage);

        //?????????????????????????????
        double min_val, max_val;
        int minLoc[2] = {255, 255}, maxLoc[2] = {255, 255};
        minMaxIdx(filtedImage, &min_val, &max_val, minLoc, maxLoc);//?????????¦Ë??maxLoc


        //????????
        Point2f center(maxLoc[1], maxLoc[0]);

        putText(greenChannel, "+", center, 0, 0.5, Scalar(0, 0, 0), 1);
        //imshow("??Center??",greenChannel);

        //????????ROI
        Mat greenROI, redROI;
        greenROI = greenChannel(
                Rect(
                        max(center.x - 80, 0),
                        max(center.y - 80, 0),
                        min(center.x + 80,greenChannel.cols-1)-max(center.x - 80, 0),
                        min(center.y + 80,greenChannel.rows-1)-max(center.y - 80, 0)));
        redROI = redChannel(
                Rect(
                        max(center.x - 80, 0),
                        max(center.y - 80, 0),
                        min(center.x + 80,greenChannel.cols-1)-max(center.x - 80, 0),
                        min(center.y + 80,greenChannel.rows-1)-max(center.y - 80, 0)));
        //imshow("redROI",redROI);

//        LOG(INFO) << '2';


        //????????????????
        //
        //
        //greenROI=greenROI*1.1;
        Mat equalize, mask;
        equalizeHist(greenROI, equalize);
        //imshow("equalize",equalize);
        threshold(equalize, mask, 20, 255, 1);
        //imshow("mask",mask);
        Mat element3 = getStructuringElement(MORPH_ELLIPSE, Size(3, 3));
        dilate(mask, mask, element3);//??mask???§Ú??
        //imshow("openMask",mask);
        Mat inpaintImage;
        inpaint(redROI, mask, inpaintImage, 1, INPAINT_NS);
        //imshow("inpaintImage",inpaintImage);

//        LOG(INFO) << '2';


        //???????????????????
        Mat element1 = getStructuringElement(MORPH_ELLIPSE, Size(10, 10));//????????
        Mat element2 = getStructuringElement(MORPH_ELLIPSE, Size(5, 5));

        Mat ROI2, ROI3;

        morphologyEx(inpaintImage, ROI2, MORPH_CLOSE, element1);//?????
        morphologyEx(ROI2, ROI3, MORPH_OPEN, element2);//??????
//        LOG(INFO) << '2';

        vector<Point> seed;
        int gradient, temp;
        for (int i = 0; i < 20; i++) {
            seed.push_back(Point(80, 80));
            gradient = 0;
            for (int j = 0; j <
                            70; j++) {
                Point p1 = pol2cart(i, j);
                Point p2 = pol2cart(i, j + 1);

                if (ROI3.at<uchar>(p2) < 10)//??p2???????§³??10????????????mask?????§Ò??????
                {
                    temp = 0;
                } else {
                    temp = ROI3.at<uchar>(p1) - ROI3.at<uchar>(p2);//????????????????????
                }

                if (temp > gradient) {
                    gradient = temp;
                    seed[i] = p1;
                }
            }
            putText(ROI3, "x", seed[i], 0, 0.5, Scalar(0, 0, 0), 1);//????§Ò???seed
        }
//        imshow("????????????", ROI3);

        //??????
        vector<double> v;

        RotatedRect box = fitEllipse(seed);

        //??box?????????ROI?§Ö??80??80????????????????
        box.center = box.center - Point2f(80, 80) + center;
        //cout << box.center <<  endl << box.size << endl;

//        ofstream SaveFile("optic.txt");
//        SaveFile << "?????????????\n";
        Mat grey;
        cvtColor(img, grey, CV_BGR2GRAY);
        int count = 0;//?????§¹??????????
        int optic_count = 0;//?????????????????
        float a = box.size.height / 2+5;
        float b = box.size.width / 2+5;
        float center_x = box.center.y;
        float center_y = box.center.x;
//        SaveFile << "?????????????" << box.center << "\n";
//        SaveFile << "??????????????????\n";
//        LOG(INFO) << '2';
        for (int y = 0; y < grey.rows; y++) {
            for (int x = 0; x < grey.cols; x++) {
                if (((y - center_x) * (y - center_x)) / (a * a) + (x - center_y) * (x - center_y) / (b * b) <= 1) {
//                    SaveFile << "(" << y << "," << x << ")\n";
                    optic_count++;
                    if (grey.at<uchar>(y, x) > 10)//10?????????????????????????¦¶????????????
                    {
                        grey.at<uchar>(y, x) = 255;
                        count++;//?????§¹??????????????
                    }
                }
            }
        }
        return grey;
    }


////////////////////////////////////////////////////////////////////////////

#if (PY_VERSION_HEX >= 0x03000000)

    static void *init_ar() {
#else
        static void init_ar(){
#endif
        Py_Initialize();

        import_array();
        return NUMPY_IMPORT_ARRAY_RETVAL;
    }

    BOOST_PYTHON_MODULE (pbcvt) {
        //using namespace XM;
        init_ar();

        //initialize converters
        to_python_converter<cv::Mat,
                pbcvt::matToNDArrayBoostConverter>();
        pbcvt::matFromNDArrayBoostConverter();

        //expose module-level functions
//        def("dot2", dot);
//        def("dot2", dot2);
        def("vesselseg", vessel_seg);
        def("findPeak", findPeak);
        def("exudate_extract", exudate_extract);
        def("hemorrhage_extract", hemorrhage_extract);
        def("getOpticDisc", getOpticDisc);

    }

} //end namespace pbcvt
